# -------------------------------------------------------------------
# Purpose: Creates Figure B12
# Author:  Max Posch, 25/07/2025
# Usage:   Source this script to generate both panels of the figure
# -------------------------------------------------------------------
# Check that required paths exist
stopifnot(dir.exists(pdataanalysis))
stopifnot(dir.exists(poutputappendix))


## Load data
load(file.path(pdataanalysis, "countyLevel19001940.RData"))

# Left panel: Patents per 1,000 people
temp <- copy(countyLevel19001940) %>% drop_na(entropy_namelast_mp_adjp_ws, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws, iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws)
temp[, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws := resid(feols(iv_lo_entropy_namelast_mp_adjp_fe_immig_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | gisjoin_1900 + statefip^year, data = temp))]
temp[, entropy_namelast_mp_adjp_ws := resid(feols(entropy_namelast_mp_adjp_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | gisjoin_1900 + statefip^year, data = temp))]
est1 <- binsreg(entropy_namelast_mp_adjp_ws, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws, data = temp)
result <- est1$data.plot$`Group Full Sample`
p <- ggplot() +
    labs(
        x = "Predicted surname diversity",
        y = NULL,
        subtitle = "Surname diversity",
    ) +
    geom_point(data = result$data.dots, aes(x = x, y = fit), size = 2, alpha = 1) +
    theme_minimal()
plotname <- file.path(poutputappendix, "figureB12left.pdf")
ggsave(plotname, p, width = 4.6, height = 4.6, units = "in")

cat("Figure saved to:", plotname, "\n")


# Right panel: Breakthrough patents per 1,000 people
temp <- copy(countyLevel19001940) %>% drop_na(entropy_namelast_mp_adjp_ws, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws, iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws)
temp[, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws := resid(feols(iv_lo_entropy_namelast_mp_adjp_fe_immig_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | gisjoin_1900 + statefip^year + gisjoin_1900[year], data = temp))]
temp[, entropy_namelast_mp_adjp_ws := resid(feols(entropy_namelast_mp_adjp_ws ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_ws | gisjoin_1900 + statefip^year + gisjoin_1900[year], data = temp))]
est1 <- binsreg(entropy_namelast_mp_adjp_ws, iv_lo_entropy_namelast_mp_adjp_fe_immig_ws, data = temp)
result <- est1$data.plot$`Group Full Sample`
p <- ggplot() +
    labs(
        x = "Predicted surname diversity",
        y = NULL,
        subtitle = "Surname diversity",
    ) +
    geom_point(data = result$data.dots, aes(x = x, y = fit), size = 2, alpha = 1) +
    theme_minimal()
plotname <- file.path(poutputappendix, "figureB12right.pdf")
ggsave(plotname, p, width = 4.6, height = 4.6, units = "in")

cat("Figure saved to:", plotname, "\n")
